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Abstract 

A  three-dimensional  mathematical  model  for  a  planar  SOFC  was  constructed.  The  concentrations  of  the  chemical  species,  the  temperature 
distribution,  the  potential  distribution,  and  the  current  density  were  calculated  using  a  single-unit  model  with  double  channels  of  co-flow  or 
counter-flow  pattern.  The  finite  volume  method  was  employed  for  the  calculation,  which  is  based  on  the  fundamental  conservation  laws  of 
mass,  energy,  and  electrical  charge.  The  internal  or  external  steam-reforming,  the  water-shift  reaction,  and  the  diffusion  of  gases  in  the 
porous  electrodes  were  taken  into  the  model.  The  effects  of  the  cell  size,  the  operating  voltage  and  the  thermal  conductivity  of  the  cell 
components  on  the  calculated  results  were  investigated.  From  the  simulated  temperature  distributions  in  the  electrolyte  and  the  inter¬ 
connector,  the  stress  distributions  were  calculated  using  the  finite  element  method.  The  results  demonstrated  that  the  steam  reforming  would 
generate  internal  stresses  of  several  tens  MPa  in  an  electrolyte.  ©  2001  Elsevier  Science  B.V.  All  rights  reserved. 
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1.  Introduction 

The  solid  oxide  fuel  cells  (SOFCs)  are  promising  candi¬ 
dates  for  future  energy  conversion  systems  because  of  their 
higher  energy  conversion  efficiency  than  those  for  conven¬ 
tional  heat  engine  systems  and  other  types  of  fuel  cells. 
In  addition,  internal  reforming  of  natural  gas  can  be  per¬ 
formed  at  the  anode,  which  allows  direct  use  of  natural  gas 
as  the  fuel. 

In  SOFC  stack  designs,  the  planar  type  design  has 
received  much  attention  recently,  because  it  is  simpler  to 
fabricate  and  easier  to  be  made  into  various  shapes  than  the 
other  type  designs.  Besides,  the  planar  type  SOFC  offers 
higher  power  density  relative  to  the  tubular  type  SOFC, 
which  is  ascribed  to  the  low  electrical  resistance  due  to 
shorter  current  paths.  However,  the  status  of  the  develop¬ 
ment  is  still  at  a  module  level,  and  some  technical  issues 
have  not  been  solved  yet.  The  internal  stresses  in  cell 
components  arising  from  thermal  shocks  or  heat  cycles 
are  one  of  the  problems  to  be  solved.  The  planar  type  SOFC 
requires  high  temperature  gas  seals  at  the  edges  or  around 
the  internal  gas  manifolds.  For  this  purpose,  cement,  glass 
and  glass-ceramic  seals  are  expected  to  give  sufficient 
sealing  efficiency.  However,  the  strict  binding  among  each 
cell  component  generates  mechanical  constraints,  and  thus  a 
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slight  mismatch  in  thermal  expansion  coefficient  (TEC) 
among  the  cell  components  can  cause  a  large  stress.  More¬ 
over,  because  of  non-homogeneous  temperature  distribu¬ 
tions  inside  the  cell,  cell  components  are  irregularly 
deformed,  adding  a  large  internal  stress.  As  a  consequence, 
the  thermal  expansion  coefficient  matching  among  the  cell 
components  and  mitigation  of  the  non-homogeneous  tem¬ 
perature  distributions  in  the  cell  components  are  indispen¬ 
sable  to  reduce  the  internal  stress.  Some  parameters  of  stack 
design  and  operating  conditions,  such  as  cell  geometry, 
internal  or  external  steam  reforming,  thermal  conductivity 
of  the  components,  and  cell  operating  voltage  affect  the 
temperature  distributions  inside  the  cells  in  a  complicated 
way,  and  thus  it  is  difficult  to  optimize  these  parameters 
independently  to  lower  the  internal  stresses.  The  computer 
simulation  technique  has  been  used  to  investigate  the  cell 
performance  in  SOFCs,  and  considered  to  be  an  effective 
way  to  analyze  the  effects  of  these  parameters  on  the 
temperature  distributions  inside  the  cells  [1-8].  However, 
only  a  few  analyses  on  the  internal  stresses  have  been 
developed  from  the  cell  performance  analyses  [9,10] 
because  fully  three-dimensional  (3-D)  model  is  required 
for  the  3-D  stress  calculation. 

The  purpose  of  the  present  work  is  to  estimate  the  thermal 
stresses  in  the  cell  components  using  a  3-D  model  simula¬ 
tion.  In  order  to  derive  the  stresses,  we  calculated  the 
temperature  distributions  in  the  cell  components  in  the 
first  step.  We  employed  simplified  single-unit  models  with 
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Nomenclature 

Fhj  diffusional  energy  flux  in  direction  Xj 

~g  determinant  of  metric  tensor 

h  ht  +  mmHm 

ht  thermal  enthalpy 

HA  total  enthalpy  of  air  in  the  channel 

HF  total  enthalpy  of  fuel  in  the  channel 
Hm  heat  of  formation  of  constituent  m 

Hs  enthalpy  in  the  electrolyte 

m„,  mass  fraction  of  mixture  constituent  m 
Ma  total  mass  of  air  in  the  channel 
MF  total  mass  of  fuel  in  the  channel 
n  unit  vector  normal  to  the  boundary 

p  piezometric  pressure 

q  heat  loss  term 

Q  heat  source  term 

,vh  energy  source 

S-,  momentum  source  components 

S{  rate  of  production  or  consumption  of  species  i 
t  time 

T  temperature 

Uj  absolute  fluid  velocity  component  in  direction  x,- 
Uj  relative  velocity  between  fluid  and  local  co¬ 
ordinate  frame  that  moves  with  velocity  ucj 
Xj  Cartesian  coordinate  ii  —  1,  2,  3) 

Greek  letters 

X  thermal  conductivity 

k  heat  transfer  coefficient 

p  density 

x jj  stress  tensor  component 


co-flow  or  counter-flow  bipolar  channels,  and  calculated  the 
concentrations  of  the  chemical  species,  the  potential  dis¬ 
tribution,  the  local  and  average  current  density,  and  the 
temperature  distribution  inside  the  cell.  The  simulation 
model  was  constructed  in  detail  to  include  the  internal  or 
external  steam-reforming,  the  water-shift  reaction,  and  the 
diffusion  of  the  gaseous  species  in  the  porous  electrodes. 


In  the  second  step  the  tensile  stresses  in  the  cell  components 
were  computed  from  the  temperature  profiles.  The  internal 
stresses  were  estimated  as  a  function  of  the  cell  size,  the 
operating  voltage,  and  the  thermal  conductivity  of  the 
cell  components,  and  thus  suitable  operating  conditions  to 
reduce  the  internal  stresses  were  proposed. 

2.  Mathematical  model 

2.1.  Model  geometry 

In  general,  one  cell-stack  is  constructed  from  a  cell 
(anode/electrolyte/cathode)  and  an  inter-connector.  As 
guides  of  an  air  and  a  fuel  to  the  cell,  an  inter-connector 
with  channel  ribs  is  employed.  The  cell-stack  is  assembled 
using  separators  as  an  auxiliary,  and  air  and  fuel  manifolds 
are  formed  to  introduce  the  air  and  the  fuel.  Fig.  1  shows  the 
schematic  diagram  of  a  cell-stack  with  internal  manifolds 
formed  in  the  separators.  Since  a  lot  of  time  is  needed  to 
simulate  the  cell  performance  for  the  whole  stack,  a  sim¬ 
plified  single-unit  cell  model  with  bipolar  channels  was  used 
in  the  present  simulation.  Only  the  counter-flow  and  the  co¬ 
flow  model  were  considered  because  numerous  meshes  are 
necessary  to  model  the  cross-flow  pattern  and  thus  the 
calculating  time  becomes  enormous.  Further,  for  the  sake 
of  simplicity  in  the  calculation,  we  analyzed  half  of  the  one 
repeating  unit  located  in  the  middle  part  of  one  center  stack. 
The  one  cell-stack  and  the  single-unit  model  are  illustrated 
in  Fig.  2.  In  this  model,  the  thickness  of  the  electrolyte, 
cathode,  and  anode  are  0.1,  0.15,  and  0.1  mm,  respectively. 

2.2.  Thermo-fluid  model 

The  thermo-fluids  analysis  was  performed  using  the 
computational  fluid  dynamics  tool  ‘STAR-CD’  (Computa¬ 
tional  Dynamics  Ltd.).  In  STAR-CD,  the  algebraic  finite- 
volume  equations  are  solved.  The  solid  and  fluid  parts  were 
divided  into  small  discrete  meshes,  and  in  each  mesh,  the 
following  differential  equations  governing  the  conservation 
of  mass,  momentum,  and  energy  were  solved. 


Fig.  1.  The  schematic  diagram  of  the  cell-stack  with  internal  manifolds. 
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Fig.  2.  The  illustrations  of  the  one  cell-stack  and  the  single-unit  model 
used  for  the  simulation. 


•  The  mass  and  momentum  conservation  (the  ‘Navier- 
Stokes’  equations): 


1  d 

sfgdt 


(V8PUj) 


d_ 

dxj 


(pUjUj 


(1) 


•  The  enthalpy  conservation  equation: 

Id  dp  duj 

=  7is(^p)  +  "'s;+T«5t;+s‘  <2) 


The  specific  heat  of  the  fluid  was  treated  as  the  polynomial 
function  of  the  temperature,  while  the  thermal  conductivity 
of  the  fluid  and  the  solid,  and  the  specific  heat  of  the  solid 
part  were  assumed  to  be  independent  of  the  temperature. 

For  the  whole  model,  the  following  mass  and  energy 
balances  should  be  considered. 

•  The  mass  balance  of  the  total  gas: 


d  Ma 

d  t 
AM¥ 
At 


=  MT  -  M 


OUT 

A 


=  MT  -  Ml 


OUT 

F 


(3) 


•  The  enthalpy  balance: 

H?"  +  Or  ,r  -  .,r:  l 

d^=H‘»-Hr+QA~^-^T  w 

^  a+w-#t 


As  mentioned  above,  in  this  simulation  we  selected  and 
modeled  one  repeating  unit  located  in  the  middle  part  of  the 
center  stack.  At  the  middle  part  of  the  center  stack,  the 
boundaries  between  the  single-unit  and  the  adjacent  units  are 
regarded  to  be  thermally  adiabatic.  Hence,  for  the  thermal 
boundary  conditions  at  the  surfaces  of  the  solid  part  con¬ 
necting  to  the  next  units,  adiabatic  boundary  conditions  were 
employed.  In  a  practical  cell  stack,  the  fuel  and  the  air  are 
introduced  into  the  channel  through  manifolds.  For  the  cell- 
stack  illustrated  in  Fig.  1,  there  are  horizontal  guide  parts 
formed  by  the  separators  from  the  manifolds  to  the  channels. 
The  edge  parts  of  the  single-unit  faces  to  the  guide  parts  and 
the  heat  exchange  between  the  edge  parts  and  the  environ¬ 
ment  around  the  stack  is  indirect.  As  it  is  difficult  to 
determine  the  exact  boundary  condition  at  the  edge  parts, 
we  employed  the  adiabatic  boundary  condition  at  the  edge 
parts.  For  the  boundary  between  the  solid  and  the  fluid,  the 
following  continuity  condition  was  imposed: 

[A(jt)Vrs(jc)]  x  n  =  k[T{(x)  —  rs(x)]  (5) 

where  n  is  the  unit  vector  normal  to  the  boundary,  Ts  and  Tt 
being  temperature  of  the  solid  and  fluid  at  x  on  the  boundary, 
respectively. 

Cell-stacks  are  usually  operated  at  900-1000°C,  and  at 
such  high  temperatures  a  radiant  heat  exchange  would  have 
an  essential  effect  on  the  heat  exchange  inside  the  channels. 
Hence,  in  this  study,  we  investigated  the  effect  of  the  radiant 
heat  exchange  inside  the  channels  on  the  simulated  results. 
To  calculate  the  radiant  heat  exchange  inside  the  channels, 
the  boundaries  of  the  calculation  domain  were  all  subdivided 
into  contiguous  non-overlapping  areas  ‘patches’.  From  the 
nominal  center  of  each  patch  a  number  of  beams  were 
emitted  and  the  each  beam  was  traced  through  the  fluid 
until  it  intercepted  an  opposing  patch.  Assumed  grey  body 
surfaces,  for  small  areas  cL4,-  and  dA;-  shown  in  Fig.  3,  the 
amount  of  radiant  energy  transfer  to  or  from  each  patch  was 
then  calculated  from  the  radiant  transport  equation: 

Ac/j-j  =  cos  4>j  cos  4>j  1  E/Ej (Ej  -  Ej)  (6) 


Fig.  3.  The  illustration  of  the  radiant  heat  exchange  inside  the  channel. 
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where  ch/,_;  is  the  net  energy  exchange  between  patch  i  and 
patch  j,  r  is  the  distance  between  the  patches,  E,  and  Ej  are 
the  radiation  energy  of  patches  i  and  j,  respectively,  and  </>;- 
are  the  angles  between  a  normal  to  the  surface  and  the  line 
drawn  between  the  patches,  and  e,-  and  Sj  are  the  surface 
emissivity  for  the  patches  i  and  j.  Using  the  Stefan-Boltz- 
mann  law 

E  =  oT4  (7) 

where  er  is  the  Stefan-Boltzmann  constant,  Eq.  (6)  is 
rewritten  as 

dA.dA,  ,  a  ,, 

d qi-j  =  cos  </>,-  cos  0 ■ - y1  Ei£M Tt  ~Tt)  (8) 

where  7}  and  7}  are  the  surface  temperatures  at  the  patches  i 
and  j,  respectively.  For  all  the  patches  the  net  energy 
exchange  written  by  Eq.  (8)  should  be  calculated.  In  the 
present  model,  the  effect  of  the  radiant  heat  exchange  on  the 
heat  balance  inside  the  channels  is  considered  to  be  small. 
Because  the  channel  length  along  the  fluid  flow  is  extremely 
long  compared  with  the  size  of  the  channel  cross  section. 
When  dealing  with  the  radiant  heat  exchange  patches,  the 
temperatures  vary  between  high  and  low.  </>,•  and  ( pj  in  Eq.  (6) 
would  be  almost  n/2  and  thus  d q,_j  determined  by  Eq.  (6) 
would  be  a  small  value.  In  addition,  when  considering  the 
radiant  heat  exchange,  the  amount  of  time  needed  to  calcu¬ 
late  precise  result  would  equate  to  ten  times  of  that  for  the 
calculation  without  considering  the  radiant  heat  exchange. 
Hence  in  this  study  we  investigated  the  effects  of  the  radiant 
heat  exchange  on  the  simulated  temperature  distribution  in 
the  cell  only  for  one  calculating  condition,  and  in  a  para¬ 
metric  study  for  the  effects  of  the  calculating  conditions  on 
the  cell  performance,  we  neglected  the  radiant  heat 
exchanges  inside  the  channels  to  reduce  the  calculation  time. 

The  radiant  losses  from  the  unit  surface  to  the  environ¬ 
ment  were  neglected  in  this  model,  since  the  singleunit 
considered  is  far  from  the  stack  surface  and  thus  the  unit 
surface  is  regarded  to  be  thermally  adiabatic.  Accordingly 
the  whole  energy  balance  equation  in  the  steady  state  was 
simplified  as  follows: 

//out  +  Hom  =  Hm  +  Hm  (9) 

2.3.  Electrochemical  model 


The  reactions  taken  into  account  are  (at  anode)  as  follows: 
Steam  reforming  reaction: 


CH4  +  H20  <->  CO  +  3H2 

(10) 

Shift  reaction: 

H20  +  CO  <->  h2  +  co2 

(11) 

Electrochemical  oxidation: 

H2  +  02~  -►  H20  +  2e~ 

(12) 

CO  +  O2  ->  C02  +  2e“ 

(13) 

Ea  =  E,  -  RJ  j-  rfi 


n 

hi 

□ 

L 

u 

Si 

Fig.  4.  The  discrete  anode/electrolyte/cathode  unit  and  the  equivalent 
circuit. 


At  the  cathode, 

j02  +  2e~  — *  O2  (14) 

Electrochemical  reactions  were  assumed  to  be  instantaneous 
at  both  electrodes  and  occur  at  the  interfaces  between 
electrodes  and  electrolyte.  It  has  been  reported  that  the  shift 
reaction  is  fast  enough  [11]  and  in  this  model  we  assumed 
that  the  shift  reaction  is  in  chemical  equilibrium  at  any  point 
in  the  anode. 

For  the  calculation  of  the  electric  current,  we  considered  a 
discrete  anode/electrolyte/cathode  unit  shown  in  Fig.  4.  The 
direction  of  the  electric  current  path  in  the  electrolyte  was 
vertical  to  the  interfaces  between  the  electrolyte  and  elec¬ 
trodes,  and  the  in-plane  path  was  treated  as  negligible.  The 
electric  current  7,-  in  the  ;th  discrete  cell  was  determined  by 
the  operating  cell  voltage  Eo  as  follows: 

Eo  =  Ej  —  RJj  —  r\t  (15) 

where  Et  is  the  Nernst  potential,  R,  is  the  Ohmic  resistance, 
and  r]i  is  the  overpotential.  The  calculated  Nernst  potentials 
include  the  concentration  overpotential  in  itself  and  thus  /;,• 
mainly  consists  of  the  activation  overpotential.  In  some 
models  for  the  cell  performance,  the  resistance  of  the  cell 
has  been  estimated  by  considering  subdivision  of  the  cell 
into  the  cell  components  (i.e.  electrolyte,  anode,  cathode, 
and  inter-connector)  [12,13].  However,  in  a  practical  cell- 
stack,  contact  resistances  among  each  cell  component  are 
usually  not  negligible  and  the  experimentally  measured  cell 
resistance  would  contain  the  contact  resistances.  Accord¬ 
ingly,  in  the  present  model,  we  employed  the  resistance  for  a 
single-cell  (anode/electrolyte/cathode),  which  was  esti¬ 
mated  from  the  experimentally  measured  I—V  data  for  a 
single-cell.  The  measured  performance  of  the  single-cell 
modeled  in  this  simulation  is  shown  in  Fig.  5.  A  value  of 
50  mV  was  used  as  the  current-independent  term  of  and 
the  current-dependent  term  was  included  in  a  formula  of  the 
Ohmic  resistance.  Using  the  empirical  formula,  the  sheet 
resistance  r;  was  expressed  as  follows: 

n  -  -34  +  1.25  x  106  x  exp(-7 185/7)  +  °'°4  <16) 

Supposing  that  no  temperature  gradients  are  present  in  the 
discrete  anode/electrolyte/cathode  unit,  the  Nernst  potential 
is  a  function  of  the  average  temperature  7)  in  the  ;th  anode/ 
electrolyte/cathode  unit  and  the  oxygen  partial  pressure  po2 
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Current  density  (A/cm2) 

Fig.  5.  Performance  of  a  single-cell  at  800,  900,  and  1000°C. 


at  both  the  electrodes.  Accordingly  the  Nernst  potential  is 
given  as 


E,  = 


FRT \npOi(c) 
4  F  P 02(d) 


(17) 


where  R  is  the  gas  constant,  F  the  Faraday  constant,  and 
p 02(c)  and  Po^(a)  are  the  partial  pressure  of  oxygen  at  the 
electrolyte/cathode  and  electrolyte/anode  interfaces,  respec¬ 
tively. 

For  the  steam  reforming  reaction  [10]  at  the  anode,  some 
formulae  describing  the  reaction  rate  have  been  proposed 
[14,15],  We  measured  the  reaction  rate  for  the  methane 
reforming  reaction  using  a  Ni/YSZ  cermet,  and  in  the 
present  model  we  used  the  following  empirical  formula 
as  the  reaction  rate: 


r ( 1 )  =  k0  exp  Pck  PhIo  A  (18) 

where  A E  is  the  activation  energy,  kQ  a  constant,  p  partial 
pressure,  and  /).,  the  density  of  the  anode.  The  equilibrium 
constant  for  the  shift  reaction  [11]  is  expressed  as 


^ Aa-ASo  ,  AalnT  A0  Ay^2\ 

Fv  RT  R  R  2R  6R  J 

(19) 


where  AH0  is  the  heat  of  reaction,  AS  the  entropy  change,  a, 
/i,  and  y  the  coefficients  of  the  temperature-dependent  terms 
in  the  function  for  the  specific  heat. 


The  reaction  rates  of  the  shift  reaction  [11]  are  determined 
for  the  forward  and  the  backward  processes: 

r{  2)f  =  k\pcoPH2o  (20) 

r(2)b  =  kiKshiftPCo2PH2  (21) 

where  k\  is  a  constant. 

The  parameters  used  in  Eqs.  (18)-(21)  are  listed  in 
Table  1. 

The  electric  current  density  J  is  expressed  by  Faraday’s 
law: 


dCb  df 

J  =  4F— -  =  2F  — 
dr  dr 


(22) 


where  d02/df  and  df/dr  are  the  molar  flux  rates  of  oxygen 
and  fuel  at  the  cathode  and  the  anode,  respectively. 
Matsuzaki  and  Yasuda  have  reported  that  the  electroche¬ 
mical  oxidation  rate  of  H2  is  about  two  times  higher  than 
that  of  CO  [16].  As  a  consequence,  in  this  model,  we 
assumed  that  FI2  consumption  rate  was  double  of  CO 
consumption  rate.  If  the  shift  reaction  is  fast  enough, 
the  shift  reaction  will  be  in  the  equilibrium  state  in  the 
anode,  and  as  a  result,  the  assumption  for  the  ratio  of  the 
H2  consumption  to  the  CO  consumption  will  not  affect 
the  simulation. 

The  diffusion  of  the  gaseous  species  in  the  porous  elec¬ 
trodes  was  considered  by  specifying  the  molecular  diffu- 
sivity  for  the  each  specimen  and  specifying  the  momentum 
source  term  in  a  linearized  form.  The  detail  of  the  calcula¬ 
tion  for  the  diffusion  of  the  gaseous  species  in  the  porous 
electrodes  was  reported  elsewhere  [17]. 


2.4.  Thermal  stress  calculation 


The  principal  stress  in  the  solid  part  was  calculated  from 
the  simulated  temperature  distribution  in  the  solid  part  using 
the  finite  element  program  ‘ABAQUS’  (Hibbit,  Karlsson  and 
Sorensen  Inc.).  In  the  stress  calculation,  the  discrete  3-D 
meshes  used  for  the  thermo-fluids  analysis  were  employed, 
and  the  element  type  was  3-D  solid  elements  of  eight-node 
linear  brick.  The  calculated  temperatures  were  read  into  the 
ABAQUS  at  the  nodes  and  they  were  then  interpolated  to  the 
calculation  points  within  elements.  In  the  stress  calculation, 
it  was  assumed  that  no  mechanical  loads  were  put  on  the 
cell  components  and  the  cell  components  could  deform 
freely,  i.e.  no  constraint  condition.  The  cell  operating 
conditions  and  parameters  used  for  the  simulation  are  listed 
in  Table  2. 


Table  1 

Parameters  for  the  chemical  reactions 


Da  (kg/m3) 

k0  (mol/kg) 

A E  (J/K  mol) 

AH0  (J) 

AS  (J) 

Act  (J) 

Afi  (J) 

Ay  (J) 

k\  (mol/m3  s) 

9.3  x  102 

1.09  X  1010 

1.91  x  105 

-4.179  x  104 

-4.337  x  101 

-9.306  x  10~‘ 

2.382  x  10~2 

-1.220  x  10~5 

1.2  x  104 
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Table  2 


The  cell  operating  conditions  and  parameters  used  for  the  simulation 


Sample 

number 

Fuel 

Air 

Channel 
length  (mm) 

Thermal  conductivity  (W/m/K) 

Flow  pattern 

Operating 
voltage  (V) 

Velocity 

(m/s) 

Temperature 

(K) 

Composition3 

Velocity 

(m/s) 

Temperature 

(K) 

Separator 

Electrolyte 

i 

0.4 

1273 

2 

3 

1123 

100 

1.8 

2.25 

counter-flow 

0.7 

2 

0.57 

1273 

Pre-reformed 

3 

1123 

100 

1.8 

2.25 

counter-flow 

0.7 

3 

0.4 

1273 

2 

3 

1123 

too 

10 

2.25 

counter-flow 

0.7 

4 

0.4 

1273 

2 

3 

1173 

100 

1.8 

2.25 

co-flow 

0.7 

5 

1.2 

1273 

2 

9 

1123 

300 

1.8 

2.25 

counter-flow 

0.7 

6 

0.4 

1273 

2 

3 

1123 

100 

1.8 

2.25 

counter-flow 

0.65 

aThe  ratio  H20/CH4. 


3.  Results  and  discussion 

3.1.  Thermo-fluid  analysis 

We  first  specified  the  standard  cell  operating  conditions 
and  parameters  for  the  counter-flow  case  so  that  the  fuel 
utilization  would  be  around  80%  (the  flow  rate  of  air  was  set 
so  that  the  oxygen  utilization  would  be  around  30%).  With 
the  selected  parameters  for  the  operating  conditions,  the 
concentrations  of  the  gaseous  species,  the  Nernst  potential, 
the  distribution  of  the  current  density,  the  average  current 
density,  the  overall  fuel  utilization,  and  the  temperature 
distribution  were  calculated.  Fig.  6(a)-(d)  show  the  calcu¬ 
lated  concentration  profile  for  CH4,  H2,  CO,  and  H20  at  the 
anode,  respectively.  Here  it  should  be  noted  that  in  this 
calculation  we  neglected  the  radiant  heat  exchange  inside 
the  channels.  It  can  be  seen  that  most  of  the  supplied  CH4 


is  reformed  within  the  first  10  mm  from  the  fuel  inlet, 
suggesting  that  the  steam  reforming  reaction  is  much  faster 
than  the  velocity  of  the  fuel  flow.  Consequently,  H2  and  CO 
concentrations  rapidly  increase  within  this  area,  and  then 
gradually  decrease  as  a  result  of  consumption  by  the  elec¬ 
trochemical  oxidation  while  flowing  towards  the  fuel  outlet. 
In  the  area  underneath  the  inter-connector  rib,  the  gaseous 
species  are  transported  only  by  diffusion  in  the  porous 
electrodes.  This  leads  to  the  lower  H2  and  CO  and  higher 
H20  concentrations  than  in  the  channel  area.  Although  in  the 
present  analyses  we  did  not  study  the  effect  of  the  rib  width 
on  the  cell  performance,  the  increase  of  the  overall  potential 
drop  with  the  increase  of  the  rib  width  has  been  reported 
[18].  At  the  fuel  outlet,  80%  of  the  supplied  fuel  has  been 
consumed. 

At  just  a  short  distance  from  the  fuel  inlet,  the  cell 
temperature  drops  rapidly  because  of  the  endothermic  steam 


Fig.  6.  The  calculated  profiles  of  CH4,  H2,  CO,  and  H20  concentrations  in  the  anode  for  the  standard  case  (no.  1). 
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Temperature 


at  electrolyte 


Fig.  7.  The  simulated  temperature  profiles  in  the  cell  and  at  the  electrolyte 


for  the  standard  case  (no.  1)  without  a  consideration  of  the  radiant  heat  exchange. 


reforming  reaction  as  shown  in  Fig.  7.  In  the  middle  area 
where  the  steam  reforming  reaction  is  complete,  the  cell 
temperature  rises  because  of  the  exothermic  electrochemical 
reaction.  At  the  downstream  of  the  fuel,  the  cell  temperature 
decreases  along  the  fuel  stream,  since  the  heat  generation  is 
small  and  in  addition,  the  hot  fuel  is  cooled  by  the  incoming 
air.  Here  we  comment  on  the  effect  of  the  radiant  heat 
exchange  inside  the  channels  on  the  simulated  temperature 
distribution  in  the  cell.  Fig.  8  shows  the  temperature  dis¬ 
tribution  at  the  electrolyte  calculated  with  considering  the 
radiant  heat  exchange  inside  the  channels.  As  comparing 
Fig.  8  with  Fig.  7,  we  can  see  two  characteristic  features. 
One  is  that  when  the  radiant  heat  exchange  is  taken  into 
consideration,  the  temperature  distribution  in  the  cell 
becomes  flat.  This  comes  from  that  the  radiant  heat 
exchange  occurs  when  the  temperatures  in  the  channels 


differ  between  high  and  low  temperature  volumes.  As 
compared  the  maximum  temperature  in  the  cell  for  s  =  0 
and  e  =  0.8,  the  difference  in  the  temperature  is  30  K.  The 
value  of  30  K  equates  about  10%  of  the  difference  in  the 
temperature  along  the  fuel  flow  at  the  electrolyte.  A  further 
difference  is  that  when  the  radiant  heat  exchange  is  con¬ 
sidered,  the  area  where  the  temperature  is  the  maximum 
moves  to  the  fuel  downstream.  When  the  radiant  heat 
exchange  is  not  considered,  the  areas  with  the  maximum 
and  the  minimum  temperature  are  very  close  as  we  can  see  in 
Fig.  7.  When  the  radiant  heat  exchange  is  considered,  a 
radiant  heat  loss  at  the  hottest  area  is  larger  because  the 
temperature  is  higher  and  the  radiation  angle  (f>  is  close  to 
zero  as  shown  in  Fig.  9.  With  increasing  the  distance  from 
the  fuel  inlet  where  the  temperature  is  the  minimum  to  the 
area  with  the  maximum  temperature,  the  value  of  the  radiant 


Temperature  [  K  1 


1455 

1424 

1392 

1361 

1330 

1299 

1267 

1236 

1205 

1173 

1142 


Fig.  8.  The  simulated  temperature  profiles  in  the  cell  and  at  the  electrolyte  for  the  standard  case  (no.  1)  with  a  consideration  of  the  radiant  heat  exchange  of 
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Fig.  9.  The  illustration  of  the  radiant  heat  exchange  between  areas  with 
high  and  low  temperatures  near  the  fuel  inlet. 


(a)  Potential  [  V  ]  (b)  Current  Density  [  A/cm2  J 

Fig.  10.  The  calculated  Nemst  potential  profile  (a)  and  current  density 
profile;  (b)  at  the  electrolyte/anode  interface  for  the  standard  case  (no.  1). 

heat  exchange  becomes  small.  Thus  the  area  with  the 
maximum  temperature  may  shit  to  the  downstream  of  the 
fuel  flow  when  the  radiant  heat  exchange  is  considered. 

Fig.  10(a)  and  (b)  indicate  the  Nernst  potential  and  current 
density  profile  at  the  electrolyte/anode  interface,  respec¬ 
tively.  The  Nernst  potential  decreases  along  the  fuel  stream, 
reflecting  the  decrease  of  the  H2  and  CO  concentrations. 
Reflecting  the  lack  of  H2  and  CO,  the  Nernst  potential  is 
lower  in  the  area  underneath  the  inter-connector  rib  than  in 


the  channel.  At  the  fuel  inlet,  instead  of  the  high  potential, 
the  current  density  is  not  so  high  because  the  cell  tempera¬ 
ture  is  low  and  thus  the  cell  resistance  is  high.  Along  the  fuel 
flow,  the  current  density  once  increases  and  reaches  a 
maximum  value  of  about  0.9  A/cm2.  At  the  fuel  down¬ 
stream,  both  the  cell  temperature  and  Nernst  potential  is 
lower,  accordingly  the  current  density  become  a  lower 
value.  Although  the  maximum  current  density  is  as  high 
as  0.9  A/cm2,  the  minimum  current  density  value  is  as  low  as 
0.073  A/cm2  in  the  area  underneath  the  inter-connector 
rib  and  thus  the  average  current  density  comes  down  to 
0.39  A/cm2. 

Next  we  focus  on  the  effects  of  the  gas  flow  pattern.  When 
the  air  inlet  temperature  in  the  co-flow  pattern  is  the  same  as 
in  the  standard  case  (counter-flow  pattern),  the  drop  of  the 
electrolyte  temperature  due  to  the  endothermic  steam 
reforming  is  larger,  resulting  in  an  increase  of  the  sheet 
resistance.  Accordingly,  the  fuel  utilization  becomes  smaller 
than  in  the  standard  case.  In  order  to  increase  the  fuel 
utilization  to  that  in  the  standard  case,  the  air  inlet  tempera¬ 
ture  has  been  modified  to  be  1173  K.  All  other  operating 
conditions  are  the  same  as  in  the  standard  case.  Fig.  11 
exhibits  the  temperature  profile  in  the  cell  and  at  the 
electrolyte  for  the  co-flow  case.  In  contrast  to  the  standard 
counter-flow  case,  the  electrolyte  temperature  increases 
monotonically  from  the  fuel  inlet  to  the  fuel  outlet.  As  a 
result,  the  potential  keeps  a  high  value  from  the  fuel  inlet  to 
the  middle  area  of  the  fuel  flow,  and  decreases  rapidly  near 
the  fuel  outlet  as  shown  in  Fig.  12(a).  Although  the  current 
density  profile  shows  a  maximum  value  in  the  middle  area 
as  well  as  in  the  standard  case,  the  difference  between  the 
maximum  and  the  minimum  current  densities,  0.46  A/cm2, 
is  smaller  than  0.88  A/cm2  in  the  standard  case  as  can  be 
seen  in  Fig.  12(b). 

All  the  simulated  results  for  the  selected  cases  in  Table  2 
are  summarized  in  Table  3.  To  validate  the  numerical 
simulation  for  the  cell  performance,  a  benchmark  test  on 


Temperature  [  K 

1352 
1323 


Fig.  11.  The  simulated  temperature  profile  in  the  cell  and  at  the  electrolyte  for  the  co-flow  case  (no.  4). 
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(a)  Potential  l  V  |  (b)  Current  Density  l  A/cm2  ] 

Fig.  12.  The  calculated  Nernst  potential  profile  (a)  and  current  density 
profile;  (b)  at  the  electrolyte/anode  interface  for  the  co-flow  case  (no.  4). 


the  numerical  model  for  the  cell  performance  has  been 
carried  out  by  International  Energy  Agency  (IEA)  partici¬ 
pants  [12].  Since  the  calculating  conditions  (e.g.  the  tem¬ 
peratures  of  the  introduced  fuel  and  air,  the  fuel 
composition,  the  single-cell  properties,  the  cell  operating 
voltage  and  thus  the  fuel  utilization)  in  the  present  simula¬ 
tion  are  different  from  those  in  the  IEA  benchmark  test,  the 
present  results  cannot  be  compared  quantitatively  with  the 
results  obtained  by  IEA.  However,  the  present  results  such  as 
the  concentration  profile  for  the  gaseous  species,  and  the 
reduction  of  the  temperature  difference  in  the  cell  by 
employing  the  co-flow  pattern  are  qualitatively  agree  with 
the  previously  reported  results  [12], 

3.2.  Stress  calculation 

We  calculated  the  thermally  induced  stresses  in  the 
electrolyte  and  the  inter-connector  for  the  six  cases  in 
Table  2.  The  distributions  of  the  principal  stresses  in  the 
electrolyte  and  the  inter-connector  were  calculated  from  the 
simulated  temperature  distributions.  The  results  for  the 
standard  case  are  shown  as  an  example  in  Fig.  13.  We 
can  see  a  stress  concentration  near  the  fuel  inlet  coming 
from  the  steep  temperature  drop  caused  by  the  endothermic 
steam  reforming  reaction.  The  maximum  principal  stress  in 
the  inter-connector  is  about  2  MPa,  which  is  low  enough  to 
avoid  damage  in  the  inter-connector.  On  the  other  hand, 
the  principal  stress  in  the  electrolyte  is  almost  70  MPa. 


Fig.  13.  The  tensile  stress  profiles  in  the  electrolyte  (a)  and  the  inter¬ 
connector;  (b)  for  the  standard  case  (no.  1). 


When  8YSZ  is  used  as  the  electrolyte,  the  internal  stress 
would  cause  cracks  or  destruction  of  the  electrolytes. 
Because  in  our  measurement  for  the  tensile  strength  on 
the  tape-cast  8YSZ,  the  compensated  tensile  strength  of 
the  electrolyte  sheets  at  room  temperature  were  estimated  to 
be  about  80  MPa  [19],  and  it  has  been  reported  that  the 
mechanical  strength  of  8YSZ  becomes  low  at  high  tempera¬ 
tures  [20,21],  All  the  calculated  maximum  tensile  stresses  in 
the  electrolyte  and  the  inter-connector  for  the  six  cases  listed 
in  Table  2  are  summarized  also  in  Table  3.  For  all  the  cases. 


Table  3 

Comparison  of  the  simulated  results 


Sample 

number 

Fuel 

utilization  (%) 

Average  current 
density  (A/cm2) 

Cell  temperature  (K) 

Maximum  tensile  stress  (MPa) 

At  fuel  outlet 

At  air  outlet 

Maximum 

In  separator 

In  electrolyte 

i 

79 

0.398 

1213 

1326 

1488 

1.96 

67 

2 

77 

0.308 

1200 

1511 

1536 

0.57 

7.7 

3 

82 

0.41 

1291 

1320 

1398 

0.45 

79 

4 

77 

0.386 

1344 

1338 

1352 

0.79 

26 

5 

77 

0.386 

1162 

1331 

1571 

2.64 

95 

6 

92 

0.463 

1232 

1443 

1617 

2.41 

77 
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Fig.  14.  The  tensile  stress  profile  in  the  electrolyte  for  the  external¬ 
reforming  case  (no.  2). 


the  maximum  tensile  stresses  in  the  electrolyte  are  several 
tens  MPa  while  those  in  the  inter-connector  are  less  than 
several  MPa. 

The  stress  concentration  would  occur  near  the  minimum 
temperature  point,  and  thus  the  magnitude  of  the  maximum 
tensile  stress  can  be  related  to  the  magnitude  of  the  steep 
temperature  drop  caused  by  the  prevailing  steam  reforming 
reaction.  For  the  external  reforming  case,  the  maximum 
tensile  stress  is  one  order  of  magnitude  smaller  than  those  for 
the  internal  methane  reforming  cases  because  there  is  no 
steep  temperature  drops  coming  from  the  largely  endother¬ 
mic  reaction.  In  addition,  the  stress  concentration  occurs  at 
the  different  region  from  that  for  the  internal  reforming  case 
as  shown  in  Fig.  14.  For  the  co-flow  case,  one  can  see  in 
Fig.  15  that  the  maximum  tensile  stress  in  the  electrolyte  is 
notably  reduced  as  compared  with  the  counter-flow  case. 
This  comes  from  the  fact  that  the  steam  reforming  reaction  is 
less  prevailing  than  in  the  counter-flow  case,  and  the  elec¬ 
trolyte  temperature  gradually  increases  along  the  fuel  flow  as 


Fig.  15.  The  tensile  stress  profile  in  the  electrolyte  for  the  co-flow  case 
(no.  4). 


Fig.  16.  The  tensile  stress  profile  in  the  electrolyte  for  the  case  of  using  an 
inter-connector  with  a  high  thermal  conductivity  (no.  3). 


shown  in  Fig.  11.  In  case  of  the  inter-connector  with  a  high 
thermal  conductivity,  despite  that  the  calculated  temperature 
difference  in  the  cell  is  smaller,  the  maximum  tensile  stress 
is  higher  than  that  for  the  standard  case  as  shown  in  Fig.  16. 
This  is  attributed  to  that  near  the  fuel  inlet  the  temperature 
gradient  along  the  direction  perpendicular  to  the  gas  flow  is 
steeper.  The  steeper  temperature  gradient  is  due  to  that  the 
highly  thermally  conducting  inter-connector  carries  the  heat 
from  the  high  temperature  region  to  the  minimum  tempera¬ 
ture  region.  Thus  these  results  indicate  that  the  steam 
reforming  reaction  is  the  substantial  origin  of  the  thermal 
stress,  and  reducing  the  induced  temperature  drop  near  the 
fuel  inlet  is  an  advantageous  way  to  avoid  the  failure  of  the 
electrolyte  sheet.  If  a  mechanical  load  exists  or  the  cell 
components  constrain  each  other,  the  maximum  stress  will 
increase  more.  It  is,  however,  possible  to  reduce  the  tensile 
stress  in  the  electrolyte  by  employing  the  co-flow  pattern. 

4.  Conclusions 

The  fluid  flow  phenomena,  including  heat  transfer,  mass 
transfer,  and  chemical  reactions  in  the  planar  SOFC  were 
simulated  for  the  simplified  single-unit  model  with  the 
bipolar  channels  for  co-flow  or  counter-flow  pattern.  The 
steam-reforming  reaction,  the  water-shift  reaction,  and  the 
diffusion  of  gases  in  the  porous  electrode  were  modeled.  The 
effects  of  the  cell  size,  the  operating  voltage,  and  the  thermal 
conductivity  of  the  cell  components  on  the  calculated  results 
were  investigated.  Using  the  finite  element  method,  the 
stress  distributions  in  the  electrolytes  and  the  inter-connec¬ 
tors  were  calculated  from  the  simulated  temperature  profiles. 
It  was  found  that  the  internal  reforming  would  induce  a  steep 
drop  of  fuel  temperature  near  the  inlet,  resulting  in  large 
tensile  stresses  in  the  electrolytes.  The  co-flow  pattern  is 
advantageous  to  mitigate  the  steep  temperature  gradient,  and 
hence  to  reduce  the  internal  stresses. 
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